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I.  INTRODUCTION 


This  report  represents  the  implementation  of  WKB  methods  [1]  for  extend¬ 
ing  ELF/VLF  modal  height  gains  in  the  ionosphere  up  to  satellite  heights.  In 
past  works  [2,3]  this  extension  has  been  acccompl ished  by  carrying  out  full- 
wave  integrations  up  to  the  satellite  altitude.  For  altitudes  sufficiently 
high  in  the  ionosphere  (ie,  >  110  km  for  VLF,  >  150  km  for  ELF  under  daytime 
ionospheres  and  i  250  km  for  ELF  under  nighttime  ionospheres)  the  electromag¬ 
netic  field  is  in  the  outgoing  whistler  mode  and  the  ionospheric  gradients  are 
sufficiently  weak  to  justify  use  of  the  WKB  method.  It  has  been  our  experi¬ 
ence  that  in  carrying  fields  up  to  altitudes  >  500  km,  the  WKB  method  will 
reduce  the  computer  cost  (Univac  1100)  by  about  a  factor  of  1/3  in  the  lower 
ELF  band  (=  75  Hz)  and  by  more  than  an  order  of  magnitude  in  the  VLF  band. 
Thus,  with  this  cost  improvement,  it  seems  quite  possible  that  case  studies  of 
satellite  reception  of  VLF  signals  using  waveguide  concepts  such  as  reported 
in  reference  4  could  be  extended  to  global  coverage  studies. 

Implementation  of  the  WKB  theory  involves  no  more  than  modifications  of  a 
full-wave  fields  program  [2]  developed  originally  for  the  purpose  of  calculat¬ 
ing  mode  conversion  coefficients  associated  with  horizontal  inhomogeneities  in 
the  earth  ionosphere  waveguide  [5].  One  departure  of  the  present  fields  pro¬ 
gram  from  that  discussed  in  reference  [2]  relates  to  the  replacement  of  the 
Inoue-Horowitz  [6]  integration  algorithm  by  the  Runge-Kutta  algorithm.  This 
change  along  with  the  capability  to  iterate  the  mode  equation  and  calculate 
excitation  factors  was  implemented  by  CH  Shellman,  of  the  NOSC  EM  Pro¬ 
pagation  Technology  Division,  in  unpublished  work.  The  Runge-Kutta  numerical 
integration  of  the  field  equations  together  with  the  Gram-Schmidt  orthogonal - 
ization  procedure  and  normalization  of  solutions  is  a  procedure  first 
described  by  Pittewqy  [7].  The  remaining  modifications  are  related  specific¬ 
ally  to  the  field  matching  required  to  Implement  the  formulas  given  In  [1], 
The  essence  of  the  method  Is  to  calculate  via  full-wave  Runge-Kutta  integra¬ 
tion  the  fields  up  to  an  altitude  called  TOPHT.  At  TOPHT  the  outgoing 
whistler  mode  (or  the  outgoing  wave  with  minimum  attenuation)  is  extracted 
from  the  set  of  four  magneto-ionic  modes  and  then  matched  to  the  corresponding 
full -wave  field  components.  The  WKB  fields  are  then  propagated  to  higher 
altitude  using  a  Simpson  rule  Integration  routine  which  adjusts  the  step  size 
in  order  to  maintain  precision. 
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Budden's  WKB  formulas  [1]  are  summarized  in  section  II.  Excitation 
factors,  useful  If  mode  sums  are  ultimately  the  goal,  and  WKB  mode  summing 
formulas  are  summarized  in  section  III.  Flow  of  the  program  and  subrouti.-e 
descriptions  are  given  In  section  IV.  Section  V  contains  a  discussion  of  card 
deck  arrangement  for  Input  along  with  discussion  of  a  sample  output.  In 
section  VI  WKB  height  gain  results  are  compared  with  full -wave  Runge-Kutta 
results,  and  the  appendix  contains  a  program  listing. 


II.  SUMMARY  OF  WKB  FORMULAS 


For  plane  wave  Incidence  of  an  rf  wave  on  the  ionosphere,  Maxwell’s 
equations  can  be  written  as  [8] 


e*  =  -iTe  (1) 

where  the  prime  denotes  (k_13/az),  with  k  the  free  space  wave  number,  T  a  4  x 
4  matrix  given  by  Budden  [1],  and  e  a  column  matrix  composed  of  components  of 
the  electric  (?)  and  magnetic  (H)  fields  of  the  rf  wave.  The  transpose  of  e 
is 


eT  =  (Ex,  -  Ey,  Z0HX,  Z0Hy)  (2) 

where  Z0  is  the  free  space  impedance.  Henceforth  the  notation  Hx  =  Zqhx,  Hy  = 
l0Hy  and  Hz  =  ZqHz  will  be  used. 

The  matrix  T  has  four  characteristic  roots  or  eigenvalues,  qi  (i  = 
1,2 ,3 ,4),  which  satisfy  the  characteristic  equation 

det  (T  -  qi)  =  0  (3) 

where  I  is  the  unit  4  x  4  matrix.  This  equation  is  the  Booker  quartic. 
Corresponding  to  any  root  q^  there  is  an  eigencolumn  P  =  s^  which  satisfies 

T  s*1*  =  qis<’>  (4) 

Let  S  be  the  4x4  matrix  whose  columns  are  s^.  For  points  where  S  is  non¬ 
singular,  the  column  matrix  f  can  be  defined  as 

e  =  Sf  or  f  =  S-^e  (5) 


and  It  can  be  shown  [8]  that  the  elements  of  f  satisfy 


f'  =  -iq  f.  +  Sr  f 
k  k  k  j  kj  J 


; j  ,k  =  1,4 


(6) 


where 
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r  =  -s'ls‘ 


(7) 


The  preceding  transformation  can  be  carried  out  only  when  S  is  non¬ 
singular.  When  two  roots  of  the  Booker  quartic  are  equal,  two  of  the  columns 
S(i)  are  usually  multiples  of  each  other,  and  then  S  is  singular.  Near  such 
points  some  of  the  coupling  coefficients  rkj  are  very  large  and  the  points  may 
be  points  of  reflection  or  points  where  coupling  between  two  upgoing  (or 
downgoing)  waves  is  very  strong.  The  present  program  cannot  be  used  in  such 
circumstances. 

When  the  species  densities  and  collision  frequencies  vary  slowly  enough 
with  height  and  where  no  two  of  the  q^  become  nearly  equal,  the  terms  of  r  are 
small  quantities  and  there  is  an  approximate  solution  for  which  the  non¬ 
diagonal  elements  of  r  are  ignored.  This  solution  is  associated  with  one 
particular  root  qj  of  the  Booker  quartic.  It  is  [I] 

f.  =  exp  (-ik/zq,dz  +  k/zr.  dz)  (B) 

J  J  J  J 

and  the  corresponding  field  components  [in  Budden's  notation]  are 

(Ex,  Ey,  Hx,  Hy)  =  ( AjFj )-l/2  +  a4,  -Aj,  qjAj,  a5qj  +  a6) 

x  exp  (-ik/zq.dz  +  k/zr,-dz)  (9) 

0  J  J 

where 

al=  '*T11  +  T44  )  a4  =  T14  T42  "  t12  t44 

a2  =  Tn  T44  -  T14  T41  a5  =  T42  (10) 

a3  =  T12  *6  =  T41  t12  "  h\  t42 

Aj  =  q-  +  Vj  +  a2 

F j  =  2qjAj  +  (q2j  -  V  Uq,  ♦  ax)  -  (a^  ♦  b^) 


(12) 

(13) 
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2rjj  ■  (Vjrl|,j  -  a3bS  +  a5b3  -  a5b3) 

+  a  b'  -  a'b  +  a  b  -  a  b  +  a  (a  b'  -a'b  +  a  b' 
4  6  4  6  6  4  6  4  j  3  6  3  6  4  5 


a.'b.  +  a_b' 
4  5  5  4 


a'b  tab'-  a'b  )1 
5  4  6  3  6  3  1 


(14) 


In  these  equations  the  T^'s  are  the  elements  of  the  T  matrix  given  by  Budden 
p  389. 

The  essence  of  the  program  documented  in  this  report  is  the  implementa¬ 
tion  of  equation  (9)  for  altitudes  exceeding  a  height  termed  TOPHT.  The  mode 
extracted  from  the  magneto- ionic  set  is  the  least  attenuated  outgoing  wave. 
Full -wave  Runge-Kutta  integrations  are  used  to  calculate  the  rf  field  compon¬ 
ents  up  to  TOPHT.  The  full -wave  solutions  at  TOPHT  are  then  matched  to  the 
corresponding  WKB  components  and  the  fields  carried  to  higher  altitudes  via 
equation  (9). 
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III.  EXCITATION  FACTORS  AND  MODE  SUMS 

Available  as  output  from  the  present  program  are  excitation  factors  for 
sub-ionospheric  electric  point  dipole  excitation  of  electric  and  magnetic  modal 
fields.  They  are  defined  as  follows: 


EXC(i.j)  =  2 


MEXC ( i ,  j )  =  2 


A4  (16) 


where 


.  s5/2 11  -  .  s3/2.ri 

1  ‘  w®  J,  1  *2 


'3  (dF/de)  ’  A4 


s‘/2  11  -  .R..R,> 

(dF/de)  R± 


In  these  equations  F  is  the  modal  function 
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and  dF/de  its  derivative  with  respect  to  the  plane  wave  angle  e.  R  and  "K  are 
plane  wave  reflection  coefficients  of  the  ionosphere  and  ground  respectively 
referenced  to  level  d.  Consistent  with  usual  notation,  the  first  subscript 
refers  to  the  polarization  (i  for  vertical  and  1  for  horizontal)  of  the  inci¬ 
dent  wave,  and  the  second  subscript  applies  to  the  polarization  of  the 
reflected  wave. 

For  the  purpose  of  defining  field  components  and  mode  sums  an  x,  y,  z 
Cartesian  coordinate  system  is  assumed  with  positive  z  directed  vertically 
upwards  into  the  ionosphere  and  x  measured  horizontally,  x-z  being  the  plane 
of  incidence.  Correspondence  between  the  i,j  indices  (see  equations  15  and 
16)  and  the  x,  y,  z  coordinate  system  is: 


i  =  j  =  1  +  z 
i  =  j  =  2  ♦  y 
i  =  j  =  3  ♦  x 


(19) 


In  terms  of  the  excitation  factors  ( EXC ( i , j ) )  the  adiabatic  (or  WKB) 

mode  sums  for  the  total  electric  field  component  Ej  may  be  written  as 

Q 

Ei(x,0,zJ= - — T7p -  S  {(EXC(l,j)T  EXC(l,j)R)J/2  cos  (Tf)ej  (z  ) 

J  K  [sin(x/a)]1/2  ln  1 

+  ( EXC ( 2 , j ) T  EXC(2,j)R>J/2  sin (7)  sin(*)  (Zj) 

+  (EXC(3,j)T  EXC(3,j)R)(|^2  sin(7)cos(0)el  (zT)) 


eR  (7  \_-ik(Sn-l)x 
ejnuR>e 


where  the  subscript  n  denotes  the  mode  index,  a  the  earth's  radius,  k  the  free 
space  wave  number  and  S  the  sine  of  the  ground  eigenangle.  The  transmitter  is 
at  (0,  0,  Zj)  and  the  receiver  is  at  (x,  0,  z^).  The  components  of  the 
electric  field  height  gain  are  denoted  by  the  ej's  and  the  superscripts  T  and 
R  signify  that  the  quantities  are  to  be  evaluated  for  the  transmitter  or 
receiver  regions  respectively.  With  the  frequency  f  in  kHz,  the  radiated 
power  P  in  kW  and 


Qe  =  6.803  x  10l  ( f kHz  PkW) 


1/2 


(21) 


,i 


the  total  field  component,  Ej,  is  in  units  of  uV/m.  The  angles  y  and  * 
measure  the  orientation  of  the  transmitting  dipole  relative  to  the  x,  y,  z 
coordinate  system  as  shown 

z 


In  similar  fashion  the  mode  sums  for  the  components  of  the  total 
magnetic  induction,  Bj,  may  be  written  as  follows 

T  R  1/2  T 

B  (X,0,ZR)  = - —  z{(MEXC(l,j)  MEXC(l,j)  )p  cosy  e^tz-j.) 

J  Csi n( x/a) ]  ' 

+  (MEXC  (2,j)TMEXC  (2,j)R)^2sin  y  sin  ^  ejn( zT)  +  (MEXC  (3,j)T 

MEXC  (3J)R)J/2  sin  y  sin  *  ejn  (zT)}hRn(zR)  e"lk^Sn_1)  x  (22) 

The  components  of  the  magnetic  field  height  gain  are  denoted  by  the  hj's. 
With 

Qb  =  7.176  x  10'2  (fkHz  Pkw)1/2  (23) 

the  total  field  component,  B  . ,  i s  units  of  y  (1y  =  10  gauss  10  W/m 

J 
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Although  the  transmitter  altitude  Zj  must  be  in  the  isotropic  portion  of 
the  guide  beneath  the  Ionosphere,  the  receiver  altitude  zR  in  equation  (20) 
or  (22)  can  be  within  the  earth- ionosphere  waveguide  or  within  the  Ionosphere 
Itself.  The  program  can,  then,  be  used  to  generate  required  inputs  for 
studies  such  as  satellite  reception  of  ground  based  VLF  transmission  using 
waveguide  concepts  [4],  Apart  from  impedance  considerations  the  program  can 
be  used  in  conjunction  with  general  magneto-ionic  reciprocity  theory  [9]  to 
estimate  fields  radiated  into  the  earth-ionosphere  waveguide  by  satellite- 
borne  sources.  However,  more  direct  means  are  available  for  making  that  type 
of  calculation  [3]  and  it  is  recommended  that  the  WKB  method  employed  here  be 
implemented  directly  into  the  program  of  reference  [3]. 

In  past  work  [2]  the  electric  and  magnetic  field  modal  height  gains  have 
been  normalized  by  equating  the  y  component  of  the  magnetic  field  to  unity  at 
the  ground.  We  will  denote  the  electric  and  magnetic  field  components  so 
normalized  by  <f.  andli..  These  quantities  are  printed  out  in  the  current 
program  under  the  listing  "HEIGHT  GAINS  WITH  HY  NORMALIZED  TO  UNITY  AT  THE 
GROUND."  The  height  gains  e ^  and  hi  entering  equations  (20)  and  (22)  are 
related  to  e,  and  It.  as  follows: 

i  i 

e!  =  (1  +  e  /e^d) 

e2  =  (1  +  RJ  e2/e2(d) 

e3  =  (1  +  e-j/e^d) 

;  barred  quantities  are  height  gain 
components  normalized  with 
^2  s  hy  =  1  at  the  9rounc1,  (24) 

hj  =  (1  +  jTrx)  T^/e^d) 

h2  =  (1  +  ,*„)  TT2  /^(d) 

h3  =  (1  +  jH:^)  Tf3  /e2(d) 

The  height  gains  e^  and  h^  are  calculated  and  printed  when  excitation 

factor  data  are  requested  (ie,  when  IEXC.NE.O).  They  are  printed  out  under 

the  label  "HEIGHT  GAINS  NORMALIZED  FOR  USE  WITH  WCB  MODE  SUMMING  FORMULAS." 

II 


It  is  these  unbarred  height  gains  which  mist  be  used  in  conjunction  with  the 
excitation  factors  defined  by  equations  (15)  and  (16)  in  the  mode  summing 
formulas  (20)  and  (22).  The  level  d  in  equation  (24)  is  the  level  where  the 
mode  equation  is  evaluated  and  is  equivalent  in  the  program  to  LWSTHT. 

IV.  PROGRAM  DESCRIPTION 

This  section  describes  the  subroutines  in  the  WKBHTG  program  listed  in 
the  appendix.  Many  of  the  subroutines  are  only  slight  modifications  of  those 
given  in  reference  [2].  The  subroutines  WKBVAR,  QGAMMA  and  DDKXMT  have  been 
added  for  the  purpose  of  implementing  the  WKB  formalism.  Principal  output  of 
the  program  is  height  gain  functions.  Height  gains  with  two  different 
normalizations  are  available.  The  first,  which  is  always  calculated,  printed 
and  plotted,  is  printed  under  the  heading  "HEIGHT  GAINS  NORMALIZED  WITH  HY 
EQUAL  TO  UNITY  AT  THE  GROUND."  This  normalization  is  consistent  with  past 
works  [2,3,5].  The  second  set  of  height  gains  is  calculated  and  printed  only 
when  excitation  factor  data  are  requested  (ie,  when  IEXC  .NE.  0).  This 
second  set  of  height  gains  is  printed  out  under  the  heading  "HEIGHT  GAINS 
NORMALIZED  FOR  USE  WITH  M<B  MODE  SUMMING  FORMULAS,"  and  defined  by  equations 
(24)  of  section  III.  Also  printed  out  when  IEXC  .NE.  0  are  the  excitation 
factors  defined  by  equations  (15),  (16)  and  (17)  of  section  III.  The  second 
set  of  height  gains  must  be  used  in  conjunction  with  the  excitation  factors 
when  implementing  the  WKB  mode  summing  formulas  (20)  and  (22).  A  chart 

showing  the  essential  structure  of  the  WKB  height  gain  program  follows  on 
pages  19  through  22. 

SUBROUTINE  MAIN 

MAIN  calls  for  the  input  of  ionic  species  data  in  XINPUT  and  for 
computation  of  height  gain  field  components  via  WAVFLD.  After  executing 
WAVFLD,  height  gains  are  available  in  the  arrays  EX,  EY,  EZ,  HX,  HY,  and  HZ  at 
DELHT  intervals  between  the  ground  and  WKBHT.  If  IEXC  .EQ.  0,  the  height 
gains  available  in  the  arrays  are  those  associated  with  HY  normalized  to  unity 
at  the  ground.  Otherwise  the  height  gains  available  in  the  arrays  are  those 
renormalized  according  to  equations  (24). 
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SUBROUTINE  XINPUT  (ISTART,  ISTOP) 

XINPUT  controls  read-in  of  input  parameters  via  NAMELIST  statements  and 
ionic  species  densities  and  collision  frequencies  as  a  function  of  altitude. 
Common  areas  are  set  up  as  required.  ISTART  is  set  to  1  before  the  first  call 
to  XINPUT  and  to  0  upon  subsequent  calls.  ISTART  =  1  Implies  all  necessary 
data  are  to  be  read  in  and  ISTART  =  0  signals  that  previously  read  data  are  to 
be  updated,  with  all  unspecified  parameters  remaining  unchanged.  If  a  value 
ISTOP  =  0  is  returned  by  XINPUT,  then  more  input  data  are  specified  in  the  data 
deck  for  subsequent  calls  to  WAVFLD,  whereas  if  a  value  ISTOP  =  1  is  returned, 
the  data  read  were  the  last  data  in  the  data  deck,  so  that  XINPUT  should  not  be 
called  again. 

The  data  deck  is  divided  into  several  parts,  each  of  which  is  marked  by  an 
identifier  card  with  the  identifier  OATUMFOL,  SPECIE,  PROFILE,  COLLFREQ,  QUIT 
or  STOP.  Each  of  these  identifiers  is  described  in  the  following  section. 

SUBROUTINE  WAVFLD  (EX,  EY,  EZ,  HX,  HY,  HZ) 

WAVFLD  calls  for  Runge-Kutta  integration  of  the  field  equations  from 
TOPHT  to  LWSTHT  at  DELHT  equally  spaced  increments  and  for  combining  the 
solutions  at  LWSTHT  so  that  they  satisfy  the  proper  modal  polarization 
condition.  It  performs  the  back  substitution  of  normalizing  values  (saved  as 
data  in  WFSTOR).  It  also  calls  for  calculation  of  height  gains  between  the 
ground  and  LWSTHT  at  DELHT  increments  via  modified  Hankel  functions  of  order 
one -third  as  well  as  the  extension  of  height  gain  functions  from  TOPHT  to 
WKBHT  via  Budden's  WKB  formalism  (these  height  gains  are  normalized  with  HY  = 

1  at  the  ground  and  are  always  printed  out).  If  IEXC.NE.0  it  also  calls  for 
calculation  of  the  excitation  factors  defined  in  section  III  and 
renormalization  of  the  height  gains  according  to  equations  (24).  Upon 
exiting,  WAVFLD  places  the  field  strengths  in  the  arrays  EX,  EY,  EZ,  HX,  HY 
HZ  at  DELHT  intervals  between  the  ground  and  WKBHT.  If  IEXC  .EQ.  0,  the 
height  gains  available  in  the  arrays  are  those  associated  with  HY  normalized 
to  unity  at  the  ground.  Otherwise  the  height  gains  in  the  arrays  are  those 
renormalized  according  to  equations  (24). 


SUBROUTINE  ITRATE 
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URATE  is  the  control  routine  for  finding  an  angle,  theta,  which 
satisfies  the  modal  equation.  It  calls  for  integration  through  the  ionosphere 
(WFINTG),  and  for  R,  1?  and  modal  function  evaluations  a  variable  number  of 
times  depending  upon  whether  iteration  of  the  input  angle  is  desired  (ie, 
ITR.NE.O)  and  whether  excitation  factors  are  desired  (ie,  IEXC.NE.O). 

SUBROUTINE  WFINTG  (TOPHT,  LWSTHT,  DELHT,  IFLAG) 

WFINTG  performs  the  Runge-Kutta  integration  of  the  two  solution  vectors 
in  P  down  through  the  ionosphere.  Numerical  solutions  are  obtained  at  all 
height  increments  of  DELHT  between  TOPHT  and  LWSTHT.  Accuracy  is  maintained 
by  continually  adjusting  the  step  size  used  in  the  numerical  integration.  The 
current  step  size  (call  it  h)  is  used  to  obtain  an  estimate  of  P,  and  then  a 
better  estimate  is  obtained  by  using  two  integrations  with  step  size  h/2.  If 
the  two  solutions  agree  within  an  error  of  PRECSN  (an  input  parameter  normally 
set  to  3.0-5),  the  better  estimate  is  accepted.  The  step  size  is  automatic¬ 
ally  decreased  to  h/Z  if  the  two  estimates  differ  by  more  than  PRECSN,  and  the 
integration  is  repeated.  If  the  error  is  significantly  greater  than  PRECSN,  a 
step  size  h/2  is  used.  If  the  error  is  significantly  less  than  PRECSN,  the 
step  size  2h  is  used  if  it  also  yields  an  error  less  than  PRECSN.  These  tests 
thus  form  an  automatic  step-size  correction.  IFLAG  is  an  internally 
controlled  flag.  If  IFLAG  =  0  the  integration  is  performed  for  THETA  only,  if 
IFLAG  =  1  the  integration  is  performed  for  THETA  and  THETA-DTHETA  where  DTHETA 
is  internally  set  to  (5.D-2,  l.D-2). 

ENTRY  INIT  T 

1NIT  T  is  an  entry  in  subroutine  TMTRX.  Sets  up  height  independent 
values  to  be  used  in  T  matrix  calculation.  These  include  the  internally  set 
flag  ISO  (ISO  =  1  for  isotropic  calculation,  0  otherwise).  Also  set  are  the 
angular  radio  frequency,  the  wave  number,  direction  cosines  of  the  geomagnetic 
field,  the  complex  sine  and  cosine  of  THETA  and  THETA-DTHETA. 

SUBROUTINE  TMTRX (HT) 

TMTRX  computes  the  value  of  the  T  matrix  at  a  specific  height  HT.  The  T 
matrix  depends  upon  input  ionospheric  parameters  (species  density , coll ision 
frequency,  angle  of  propagation,  magnetic  field,  etc).  The  susceptibility 
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matrix,  M,  for  each  species  in  the  ionosphere  is  confuted,  the  effect  .of  earth 
curvature  is  included  and  the  T  matrix  is  computed  from  the  susceptibility  matrix 
elements.  The  equations  used  to  evaluate  the  M  and  T  are  given  in  Clemmow  and 
Heading  [8]  and  by  Budden  [1]. 

SUBROUTINE  WFDENS 

WFDENS  (HT,  EN,  COLL)  computes  the  species  density  (returned  in  EN)  and 
collision  frequencies  (returned  In  COLL)  at  height  HT  for  up  to  five  species 
in  the  Ionosphere.  EN  and  COLL  are  determined  from  corresponding  profiles  in 
the  common  field  WFPROF.  LHT  and  MHT  are  integer  values  which  indicate  which 
profile  values  are  to  be  interpolated  to  find  values  at  the  height  HT.  The  EN 
(or  COLL)  values  are  given  by  logarithmic  Interpolation  of  the  profile  values 
at  heights  ENHT(LHT  +  1)  and  ENHT(LHT)  or  (COLLHT(MHT  +  1)  and  COLLHT(MHT) , ) . 

SUBROUTINE  WF  INIT(P) 

WF  INIT  computes  the  two  starting  field  vectors  in  P  at  TOPHT  subject  to 
the  condition  of  outgoing  wave.  This  is  done  by  making  use  of  Booker  quartic 
solutions  for  a  homogeneous  anisotropic  medium. 

SUBROUTINE  QUARTC  (FOUR  B3,  SIX  B2,  F0URB1,  B0,Q) 

QUARTC  computes  the  four  roots  of  the  polynomial  Q4+F0URB3*Q3  +  SIXB2*Q2 
F0URB1*Q  +  BO  In  closed  form  [10].  Up  to  ten  applications  of  Newton's 
Iteration  are  then  performed  to  improve  the  accuracy  of  the  roots,  if 
necessary. 

SUBROUTINE  PDERIV(P,DPDH) 

PDERIV  computes  the  height  derivatives  of  the  two  field  vectors  In  P 
according  to  Clenmow  and  Heading  [8].  The  derivative  Is  returned  in  DPDH. 

ENTRY  TI  MTRX 

This  Is  an  entry  In  TMTRX  which  calculates  T  matrix  elements  for  the 
incremental  plane  wave  angle  THETA-DTHETA. 


- - •  .....  ■;«!». 


SUBROUTINE  XFER  (A,B,N) 
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Transfers  the  N  element  array  A  into  B. 


ENTRY  WF  STOR 

WF  STOR  is  an  entry  in  WFSCAL  where  certain  values  obtained  during 
integration  through  the  ionosphere  are  saved  for  later  use.  The  solution 
matrix  P,  the  height  for  which  P  is  a  solution  and  a  height  integer  index  are 
saved  along  with  orthogonal izati on  and  normalization  values  OSUM,  APROD  and 
BPROD.  In  addition,  values  of  the  susceptibility  matrix  elements  M31,  M32  and 
M33,  which  are  needed  to  compute  the  EZ  component  of  the  electric  field,  are 
saved  at  each  height. 

SUBROUTINE  WF  STEP  (P.DPDH,  HT,  DELHT,  I FLAG) 

WF  STEP  numerically  advances  the  solution  matrix  P,  using  the  derivative 
DPDH,  from  HT  to  HT-DELHT  by  the  Runge-Kutta  method.  IFLAG,  set  internally, 
controls  the  calculations  at  intermediate  points  between  HT  and  HT-DELHT  at 
which  evaluations  are  required  for  comparison  of  the  second  and  fourth  order 
Runge-Kutta  integrations. 

SUBROUTINE  WF  SCAL  (PP,  IFLAG) 

WF  SCAL  normalizes  and  orthogonal izes  the  solution  vector  PP  according  to 
the  formulas  of  reference  [2].  This  scaling  must  later  be  removed  to  yield 
correct  unsealed  solutions.  Control  for  calculating  the  quantities  needed  for 
removal  of  the  scaling  is  the  internally  set  IFLAG.  The  quantities  OSUM, 
APROD  and  BPROD  needed  in  the  backward  substitution  (or  equivalently  for 
removing  the  scalings)  are  calculated  for  IFLAG  =  0. 

ENTRY  FFCT  (PP,  C,  S,  F) 

This  is  an  entry  in  WF  BNDY  which  evaluates  the  mode  function: 


F  =  (1-  R 


R  )  (1-  R  R  ) 
nit  iiii 


-iR. 


UR1  ltR  H  iRi 


The  two  field  solution  vectors  given  in  PP  along  with  cosine  (C)  and  sine 
(S)  of  the  plane  wave  incident  angle  THETA  are  used  in  computing  the 
reflection  coefficients  by  plane  wave  decomposition  at  LWSTHT. 
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SUBROUTINE  RMTRX  (P,  COSN,  R) 

RMTRX  computes  the  2x2  Ionospheric  reflection  coefficient  matrix 
(returned  In  R)  from  the  two  solution  vectors  given  in  P.  The  two  solutions 
at  LWSTHT  are  transformed  into  upgoing  and  downgoing  components  and  the  ratio 
of  these  components  determines  the  reflection  coefficients. 

SUBROUTINE  RBARS(C,S,RBAR11,  RBAR22,  EY,  HY) 

RBARS  calculates  the  plane  wave  reflection  coefficients  (returned  in 
RBAR11  and  RBAR22)  looking  towards  the  ground  from  LWSTHT.  Evaluations  are  in 
terms  of  modified  Hankel  functions  of  order  1/3.  The  cosine  (C)  and  sine  (S) 
of  the  plane  wave  angle  THETA  are  used  in  the  IT  determination.  Note  that 
RBAR11  =*  and  RBAR22  =*  T  . 

SUBROUTINE  MDHNKL  (Z,  HI,  H2,  H1PRME,  H2PRME ) 

MDHNKL  calculates  for  argument  Z  two  independent  solutions  (HI  and  H2) 
and  their  derivatives  (HIPRME,  H2PRME)  of  Stokes'  equation  by  methods 
described  in  reference  [11], 


SUBROUTINE  WF  BNDY(B) 


WF  BNDY  computes  the  coefficients  Bj,  B2  (In  B)  which  are  used  to 
linearly  combine  the  two  solutions  in  P  to  form  the  unique  solution  which 
satisfies  the  boundary  condition  on  the  modal  polarization  at  LWSTHT  and  the 
normalization  condition  HY  =  1  at  Z  =  0.  The  polarization  condition  is: 


EY/HY  = 


1  +  lRl  ^  “  lRl 


1  +  1*1 


lRl  i^i 


lRl  iRl 


I+  «R. 


1  +  ,R, 

where  all  reflection  coefficients  are  referenced  to  LWSTHT. 


1  "  iRi  1*1 


ENTRY  HT  GAIN  (ALTT,  EX,  EY,  EZ,  HX,  HY,  HZ) 

HT  GAIN  is  an  entry  in  RBARS  where  the  field  height  gains  at  heights  ALTT 
between  the  ground  and  LWSTHT  are  calculated  at  DELHT  intervals  using  modified 
Hankel  functions  of  order  1/3. 


ENTRY  WF  LOAD 


WF  LOAD  is  an  entry  in  WFSCAL  which  transfers  normalizing  and 
orthogonal izi ng  values  saved  in  WFSTOR  for  the  purpose  of  unsealing  the  two 
solutions  in  P  starting  from  LWSTHT  and  proceeding  to  TOPHT  at  DELHT 
intervals. 

SUBROUTINE  WCBVAR  (EX,  EY,  EZ,  HX,  HY,  HZ) 

WKBVAR  extends  calculation  of  the  field  equations  from  TOPHT  to  WKBHT  by 
means  of  the  formulas  in  section  II.  Integration  of  the  phase  factors  and 
T j j  (j  index  for  least  attenuated  outgoing  mode)  is  performed  by  a  Simpson 
rule  routine  which  maintains  precision  by  adjustment  of  the  step  size  much 
like  the  Runge-Kutta  step  size  adjustment  in  WFINTG. 

ENTRY  INITDT 

This  is  an  entry  in  DDKXMT  where  all  height  independent  quantities  are 
initialized  before  extending  field  calculations  from  TOPHT  to  WKBHT  via  the 
WKB  method. 

SUBROUTINE  Q  GAMMA  (HT,  DELHT,  TOPHT,  Q,  GAMMA) 

Q  GAMMA  determines  the  Booker  quartic  solution  for  the  least  attenuated 
outgoing,  magneto-ionic  component  at  TOPHT  and  computes  the  ai  1  s ,  b^s, 
a^'s,  b^’s,  A.j ,  Fi  and  according  to  formulas  of  section  II,  Full -wave 
solutions  are  matched  at  TOPHT  to  the  WKB  solutions.  At  other  heights  and  up 
to  WKBHT  coefficients  of  the  exponentials  in  equation  (9)  of  section  II  are 
calculated  along  with  coefficients  for  HZ  and  EZ. 

SUBROUTINE  DDKXMT 

DDKXMT  calculates  susceptibility  matrix  elements,  M,  the  T  matrix 
elements  and  their  derivatives  with  respect  to  height  in  the  height  range 
TOPHT  ±Z  ±  WKBHT. 

ENTRY  EXCFAC(S) 

EXCFAC  is  an  entry  in  WFBNDY  where  the  excitation  factors  for  the  sine 
(S)  of  the  plane  wave  eigenangle  THETA  are  calculated  according  to  the 
formulas  of  section  III. 
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MAIN 

Calls  for  input  data  via  XINPUT, 
computation  of  height  gain  func¬ 
tions  (with  HY=1  at  the  ground), 
excitation  factors  (if  IEXC  .NE.  0) 
along  with  renormalized  height 
gains  (if  IEXC.NE.O)  by  WAVFLD. 


XINPUT 

Input  electron  and 
ion  species  densities  and 
collision  frequencies  as 
a  function  of  altitude. 


WAVFLD 

Calls  for  Runge-Kutta  integration 
of  the  field  equations  from  TOPHT  to 
LWSTHT  at  DELHI  increments  and  for  com¬ 
bining  solutions  at  LWSTHT  so  that  they 
satisfy  the  proper  modal  polarization 
conditions.  It  performs  the  back  sub¬ 
stitution  of  normalizing  values  (saved 
as  data  in  WFSTOR).  It  also  calls  for 
calculation  of  height  gains  between  the 
ground  and  LWSTHT  at  DELHT  increments 
via  modified  Hankel  functions  of  order 
1/3  as  well  as  the  extension  of  height 
gain  functions  from  TOPHT  to  WKBHT  via 
Budden's  WKB  formalism.  These  height 
gains  are  normalized  with  HY=1  at  the 
ground.  Lastly,  if  IEXC.NE.O  WAVFLD 
calls  for  calculation  of  excitation 
factors  and  height  gains  normalized 
according  to  Eqs  (24). 
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WAVFLD 


WFINTG 


V.  SAMPLE  INPUT  AND  OUTPUT 


A.  INPUT 

Altitude  Independent  parameters,  species  densities  and  collision 
frequencies  as  a  function  of  altitude  are  supplied  via  an  input  data  deck  on  a 
standard  Input  unit.  Read -in  occurs  in  the  subroutine  XINPUT.  The  data  deck 
is  divided  into  several  parts,  each  of  which  is  marked  by  control  cards 
DATUMFOL,  SPECIE,  PROFILE,  COLLFREQ,  QUIT,  or  STOP.  Each  of  these  control 
cards  Is  described  below: 

1 )  DATUMFOL 

DATUMFOL  -  is  a  control  card  signifying  that  input  for  the 
namelist  DATUM  follows. 

&DATUM  -  initiates  reading  of  altitude  independent  input 
data  in  namelist  input  format.  These  data  are: 


THETA  -  complex  angle  of  incidence  in  degrees  as  measured 
at  H. THETA  must  be  an  eigenangle  if  ITR  =  0. 

FREQ  -  frequency  in  kHz. 

IDBG  -  Integer  flag  controlling  auxiliary  printout.  If 
IDBG.GE.l,  the  R  matrix  elements,  IT  and^.  are 
printed.  If  IDBG.GE.2  the  solution  vectors  (P)  at 
TOPHT  and  LWSTHT  are  printed  along  with  the  Booker 
quartlc  solution  at  TOPHT  and  the  combining 
coefficients  (B(l),  B{2))  of  the  two  independent 
solution  vectors  contained  in  P  at  LWSTHT. 


TOPHT  -  starting  height  for  full-wave  integration  (km). 


LWSTHT -  lowest  height  for  full -wave  integration  (km). 

Also  equivalent  to  level  D  where  mode  equation  is 
evaluated. 

WKBHT -  altitude  in  km  to  which  height  gains  are  extended 
beyond  TOPHT  via  WKB  formulas. 

DELHT  -  height  increment  in  km  at  which  height  gains  are 
printed. 

PRECSN  -  accuracy  to  be  maintained  locally  in  the  numerical 
integrations.  Usually  taken  to  be  the  default 
value  of  3.0E-5. 

AZIM  -  azimuth  of  propagation  path  in  degrees,  measured 
clockwise  from  geomagnetic  north. 

CODIP  -  codip  of  geomagnetic  field  in  degrees 

MAGFLD  _  geomagnetic  field  strength  in  webers/square  meter 
COEFNU ( 5 ) 

-coefficient  and  exponent  of  exponential  form  of 
EXPNU(5)  collision  frequency  (if  not  specified  by  profile). 

Up  to  five  values  may  be  specified,  one  for  each 
species. 

ALPHA  -  earth  curvature  coefficient  in  inverse  km. 
Default  is  3.14E-4. 

SIGMA  -  ground  conductivity  in  siemens/meter.  Default 

value  is  the  seawater  value  of  4.64. 

EPSLON  -  ground  permittivity  in  farads/meter.  Default 

value  is  7  .172015D-10.  This  corresponds  to  a 
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dielectric  constant  of  81  for  seawater. 

ISO  -  flag  signaling  whether  calculations  are  to  be 

performed  for  isotropic  or  anisotropic 
Ionosphere.  ISO.NE.O  signals  Isotropy. 

ITR  -  integer  flag  which  calls  for  modal  equation 

iteration  when  ITR  ,NE.  0. 

H  -  altitude  in  km  at  which  modified  refractive  index 

is  unity.  Elgenangles  are  referenced  to  this 
altitude. 

AEND  -  signifies  end  of  the  DATUM  namelist  input. 

11)  SPECIE 

SPECIE  -  is  a  control  card  signifying  that  input  for  namelist 
SPECIE  follows. 

&SPECIE  -  initiates  reading  of  altitude  independent  species 
data.  These  data  are: 

NRSPEC  -  number  (Integer)  of  species  in  the 

ionosphere.  Can  take  on  values  up  to  5. 
Default  value  is  1. 

CHARGE(5)  -  sign  of  charge  of  each  species  in  proton 

units.  For  an  electron,  the  CHARGE  is 
-1.0.  Default  values  are  -1.0,  1.0,  -1.0, 
1.0,  -1.0). 

RATI0M{5)  -  mass  of  each  species  relative  to  mass  of  an 

electron.  Default  values  are  (1.0,  5.8D4, 
5.8D4,  5.8D4,  5.8D4). 

AEND  -  signifies  the  end  of  the  SPECIE  namelist  Input, 

ill)  PROFILE 
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PROFILE 


is  a  control  card  initiating  reading  of  the 
ionospheric  profile  cards.  The  control  car<J 
PROFILE  is  followed  by  an  alpha-numeric  card 
which  is  used  to  identify  the  profile.  The 
profile  is  input  starting  at  the  top  of  the 
ionosphere.  The  cards  must  be  input  in 
descending  order  by  height.  The  profile 
cards  contain  the  height  in  kilometers  and 
the  species  densities  in  particles  per  cubic 
centimeter.  A  maximum  of  five  species  can  be 
specified.  In  the  special  case  of  three 
species,  only  two  are  specified  on  the 
card.  The  first  is  assumed  to  be  electrons 
and  the  second  is  positive  ions.  The  third 
species,  negative  ions,  is  calculated  by 
subtracting  the  electron  density  from  the 
positive  ion  density.  All  three  species  are 
listed  on  the  computer  printout.  If  the 
value  of  any  species  density  is  less  than  or 
equal  to  zero,  it  is  set  in  the  program  to 
1.0E-40.  The  heights  are  punched  in  the  form 
XXX .XX  with  the  decimal  point  in  column  5, 
and  the  densities  are  punched  in  the  form 
X.XXD+XX  with  the  decimal  points  in  column 
15,  25,  etc.  All  species  data  must  be 

specified  except  for  the  special  case 
discussed  above.  The  end  of  the  profile  is 
indicated  by  a  dummy  height  of  999.99. 


iv)  COLLFREQ 

COLLFREQ  -  is  a  control  card  initiating  reading  of  the 
collision  frequency  profile.  This  allows  using 
non-exponential  collision  frequencies.  A  strictly 
exponential  collision  frequency  may  be  specified  in 
namelist  input  by  the  variables  COEFNU  and  EXPNU. 
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If  a  profile  deck  Is  used  It  overrides  COEFNU  and 
EXPNU.  Collision  frequencies  for  all  species  must 
be  input.  The  card  preparation  is  just  as  above 
with  species  density  values  replaced  by  collision 
frequency  values  (in  collisions/s)  on  the  cards. 


v)  QUIT 

QUIT  -  is  a  card  control  which  indicates  the  end  of  input 

data  for  a  call  to  WAVFLD.  After  calling  WAVFLD, 
XINPUT  can  be  read  and  used  as  data  for  a 
subsequent  call  to  WAVFLD.  This  allows  several 
runs  to  be  made  with  one  input  deck.  Note  that 
only  those  data  which  are  changed  need  be  specified 
after  the  card  QUIT. 

vi)  STOP 

STOP  -  control  card  which  indicates  that  there  are  no  more 

input  data  and  that  the  run  is  to  be  tQrminated 
after  the  next  call  to  WAVFLD.  Note  that  if  QUIT 
is  encountered  I  STOP  is  set  to  zero,  but  if  STOP  is 
encountered  I  STOP  is  set  to  one. 

A  schematic  of  an  input  deck  is  shown  on  page  30  ,  and  an  actual  sample  input 

is  shown  on  pages  31  and  32. 

B.  OUTPUT 

The  output  corresponding  to  the  input  is  shown  on  pages  33  through  51. 

The  output  corresponds  to  the  flag  settings  ITR  .EQ.  1  and  IEXC  .EQ.  1. 


First  come  NAMELIST  and  profile  printouts  followed  by  output,  which  gives 
the  number  of  Runge-Kutta  integration  steps  used  in  WAVFLD  during  the  course 
of  Integrating  from  TOPHT  to  LWSTHT.  This  Is  followed  by  the  modal  equation 
values  for  the  input  eigenangle  THETA  shown  on  page  35  .  Since  ITR.EQ.l  an 
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iteration  is  performed  and  the  iterated  value  of  THETA  is  shown  as  NEW  THETA 
on  the  output.  The  process  of  integration  from  TOPHT  to  LWSTHT  is  repeated 
for  this  NEW  THETA.  Again  the  number  of  integration  steps  and  modal  equation 
value  is  shown.  The  iteration  stops  because  the  change  in  the  absolute  value 
of  the  real  part  of  the  eigenangle  is  .LE.  5.D-2  and  the  change  in  the 
absolute  value  of  the  imaginary  part  is  .LE.  5.D-3. 

The  first  set  of  field  strengths  at  50.0  km  (LWSTHT)  is  calculated  using 
ground  boundary  conditions  and  modified  Hankel  functions  of  order  one- third  to 
extend  the  calculations  to  50.  km.  The  second  set  is  calculated  from  the 
Runge-Kutta  field  integration.  Together  the  two  sets  of  field  strengths  give 
an  alternative  check  of  ho.,  well  the  modal  condition  is  satisfied.  This 
alternative  check  is  particularly  useful  when  ITR  .NE.  0  and  eigenangles  are 
input  from  another  program  such  as  the  waveguide  program  of  reference  [12]. 
It  has  been  our  experience  that  a  good  eigenangle  will  cause  the  two  sets  of 
field  values  to  agree  to  better  than  one  percent. 

Since  IEXC  .EQ.  1  the  excitation  factors  EXC(I,J)  and  MFXC ( I ,J )  defined 
in  section  III  are  calculated  and  printed.  The  rows  correspond  respectively 
to  electric  dipole  excitation  produced  by  a  vertical  (V),  horizontal  broadside 
(HB),  and  horizontal  end-on  (HE)  oriented  dipole.  The  columns  are  for 
excitation  of  the  EZ  (or  HZ),  EY  (or  HY)  and  EX  (or  HX)  components.  That  is 
to  say  EXC(I,J)  and  MEXC(I,J)  are  to  be  combined  with  height  gain  components 
according  to  the  mode  summing  formulas  of  section  III. 

Following  the  excitation  factors  come  the  height  gains  calculated  at 
OELHT  intervals  up  to  WKBHT.  They  are  normalized  with  HY  =  1  at  the  ground. 
This  is  consistent  with  past  work  [2,5].  The  magnitudes  of  these  height  gains 
are  also  plotted  and  the  plot  is  shown  on  page  52.  The  rapid  change  in  the 
slope  of  the  EY,  EZ,  and  HZ  plots  is  quite  surprising.  An  immediate  reaction 
might  be  that  this  is  where  the  full-wave  and  WKB  fields  are  matched.  How¬ 
ever,  though  unexplained,  the  break  d  94  km)  is  much  too  low  (TOPHT  =  120  km) 
and  totally  unrelated  to  field  matching  at  TOPHT.  Consistency  of  the 
full-wave  calculations,  as  indicated  by  the  modal  function  evaluation  and  the 
continuity  of  the  field  components  at  LWSTHT  along  with  the  fact  that,  as 
indicated  in  the  following  section,  the  fullwave  and  WKB  methods  yield  nearly 
Identifical  results  for  fields  above  TOPHT,  points  towards  the  accuracy  of 
the  full-wave  methods.  It  is  clear  that  very  little  of  the  EY,  EZ  and  HZ 
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components  penetrate  the  Ionosphere  for  the  first  order  mode  used  for  the 
Sample  Input.  Penetration  Is  greater  the  higher  the  order  of  the  mode. 
Finally,  the  large  value  of  the  relative  fields  within  the  guide  Is  peculiar 
to  the  low  conductivity  (10“5  S/M)  indicated  in  the  sample  input.  Because 
IEXC  .HE.  0  a  second  set  of  height  gains  normalized  according  to  Equations 
(24)  are  calculated.  These  are  the  height  gains  that  are  to  be  used  In 
conjunction  with  the  excitation  factors  EXC(I,J)  and  MEXC(I,J)  in  the  WKB  mode 
sunning  formulas  (20)  and  (22). 


EXAMPLE  OF  INPUT  DECK 


DATUMF0L 

&OATUM 
(Input  data) 

AEND 

SPECIE 

ASPECIE 

(input  specie  data) 

>END 

PROFILE 

(profile  cards  for  each  specie) 

999.99 

C01LFREQ 

(coll,  frequency  cards  for  each  specie) 

999.99 

QUIT  (end  of  input  for  this  run) 

DATUMF0L 

ADATUM 

(changes  In  input  data) 

SEND 

PR0FILE 

(specify  entire  new  profile  deck) 

999.99 

QUIT  (end  of  input  for  this  run) 

DATIJMF0L 

&DATUM 

(changes  In  Input  data) 

AENQ 

ST0P  (end  of  data  deck) 
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199  INTEGRATION  STEPS  USED  IN  WAVFLO 
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.  A 

PROGRAM  VERIFICATION 


A  variety  of  comparisons  were  made  between  the  WKB  extended  height  gains 
and  the  full -wave  Runge-Kutta  results.  Two  representative  examples  of  the 
comparisons  are  given  in  tables  1  and  2. 

Table  1  shows  comparison  between  the  full-wave  calculated  field 
components  and  the  WKB  calculations  for  the  parameters  given  in  the  sample 
input  of  the  previous  section.  Given  in  the  table  are  all  six  field 
components  at  250,350  and  500  km.  Five  sets  of  WKB  results  are  given.  These 
correspond  to  TOPHT  settings  of  100,  110,  120,  130  and  150  km.  Though  there 
is  a  reasonable  improvement  between  the  results  for  TOPHT  >_  110  Km  and  those 
for  TOPHT  =  100  Km,  the  improvement  with  increasing  TOPHT  to  150  Km  is  at  best 
very  minor.  This  is  probably  because  some  residual  reflection  occurs  at 
altitudes  greater  than  150  Km.  For  TOPHT  >_  110  Km  most  of  the  comparisons 
agree  to  better  than  2%.  There  are  a  couple  of  cases  where  the  smaller  of  the 
Re  and  Im  parts  differ  from  the  fullwave  calculation  by  about  10*.  In  view  of 
path  and  ionospheric  uncertainties  this  difference  would  be  of  no  significance 
in  any  practical  applications.  Even  with  T0PHT=150  Km,  the  cost  savings  over 
the  full  wave  integration  was  more  than  an  order  of  magnitude  on  the  Uni  vac 
1100/82  system. 
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ALTITUDE l KM)  ALTITUDE (KM)  ALTITUDE (KM) 


GLOBAL  NIGHTTIME  IONOSPHERE (SAT.  NIGHT  ABOVE  99,  h'=87  BELOW  99) 
AZIM=37 . 8  DEG, CODIP-6. 6  DEG, MAGPLD-S. 2*10  '  GAUSS 
FREO-  17.8  KhZ,  S I GMA- 1 0.  Ox  1 0'6S/M,  EPSLON-  4 .  T*  1  O'"  P/M 
AT  NORMALIZED  TO  i  AT  THE  GROUND 


MAGNITUDE  CP  EZ  MAGNITUDE  OP  hZ 


PROGRAM  VERIFICATION 


A  variety  of  comparisons  were  made  between  the  WKB  extended  height  gains 
and  the  full -wave  Runge-Kutta  results.  Two  representative  examples  of  the 
comparisons  are  given  in  tables  1  and  2. 

Table  1  shows  comparison  between  the  full-wave  calculated  field 
components  and  the  WKB  calculations  for  the  parameters  given  in  the  sample 
input  of  the  previous  section.  Given  in  the  table  are  all  six  field 
components  at  250,350  and  500  km.  Five  sets  of  WKB  results  are  given.  These 
correspond  to  TOPHT  settings  of  100,  1x0,  120,  130  and  150  km.  Though  there 
is  a  reasonable  improvement  between  the  results  for  TOPHT  ^110  Km  and  those 
for  TOPHT  =  100  Km,  the  improvement  with  increasing  TOPHT  to  150  Km  is  at  best 
very  minor.  This  is  probably  because  some  residual  reflection  occurs  at 
altitudes  greater  than  150  Km.  For  TOPHT  >_  110  Km  most  of  the  comparisons 
agree  to  better  than  2%.  There  are  a  couple  of  cases  where  the  smaller  of  the 
Re  and  Im  parts  differ  from  the  fullwave  calculation  by  about  10%.  In  view  of 
path  and  ionospheric  uncertainties  this  difference  would  be  of  no  significance 
in  any  practical  applications.  Even  with  T0PHT=150  Km,  the  cost  savings  over 
the  full  wave  integration  was  more  than  an  order  of  magnitude  on  the  Univac 
1100/82  system. 
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(See  sample  input  for  geomagnetic  and  ground  parameters.) 
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program  driver 

the  driver  program  alternately  calls 

FOR  THE  INPUT  OF  DATA  ON  THE  STANDARD 
INPUT  UNIT  AND  CALLS  FOR  THE  COMPUTATION 
OF  HEIGHT  GAIN  FUNCTIONS  BY  WAVFLD. 
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I P  V  DABS ( HT-D )  .GT. 
EX(IHT)  =  EXSTOR 
EY(IHT)  «  EYSTOR 
E2 ( IHT )  =  EZSTOR 
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CALL  MESSAG( LA8EL2. 80 ,0.8,10.1) 

CALL  MESSAG( ' AZIM=  DEG.C00IP*  DEG.MAGFLD*  GAUSS' 


IMPLICIT  REAL*8  (A-H.0-2) 

COMMON  /WFINPT/  THETA,  FREQ. 

AZMUTH,  CODIP,  MAGFLO,  CEFFNU(5),  EXPNU(5) 
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CALL  XFER(P,P  SAVE. 16) 

1 F(HT  .LT.  WFHT+EPSHT)  CALL  WF  STOR 
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